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Abstract 
We show how the error term for Gauss-Legendre quadrature can be estimated, by solving a 
suitable initial value problem using a Runge-Kutta method. The error term can then be used as 


a correction term, yielding a more accurate result for the quadrature. We have also identified a 
condition function that can amplify/suppress the effect of the Runge-Kutta global error. 


Keywords: Gauss-Legendre, Runge-Kutta, error, initial value problem, correction term 


Mathematics Subject Classification 2020: 65D30, 65L05, 65L70 


1 Introduction 


We have recently considered the possibility of determining the error term in Trapezium 
quadrature by solving a suitable initial-value problem [1]. In this paper, we extend this idea 
to Gauss-Legendre quadrature. 


2 Relevant concepts 


We establish here notation and terminology relevant to the rest of the paper. 


2.1 Gaussian Quadrature 


A Gaussian quadrature rule [2][3] has the form 


GIf.w,a,x,n+1] Cras 


i=0 
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where the parameters c;” and x,” denote the weights and nodes, respectively. A Gaussian 
rule is intended to approximate the integral 
x 
I Lf, w, a,x] = : w (x) f (x) dx, 
a 


where w(x) > 0 on [0, x] is a weight function. We necessarily assume that 
f:R-R 
w:R-R 


and that f (x) is as differentiable as it needs to be for our purposes. Naturally, we also assume 
that w (x) f (x) is integrable on any interval that we consider. 


The following statement gives the relationship between the weights, nodes and degree of 
precision of a Gaussian quadrature rule [4]: The (n + 1)-point quadrature rule 
x 


GIf,w,a,x,n+1] = =) otf (x YY oes [ re@orera 


i=0 A 


where w(x) is a Meer function on [a, x], has degree of precision 2n + 1 if and only if the 
nodes tag Re ikea et are the roots of the w-orthogonal polynomial P™ 


the weights c” are ave by 
P hf xx," 
- f wof]| — ra 
9 (Xi 3% 


a na 


W(x) on [a, x], and 


jFi 
The approximation error in G[f, w, a, x, n + 1] is given by 


x 2 
(2n+2) ( (x)) n 
f wt eae otfnaney Led f w (x) (x -xi”)] dx, 
a a = 
where C € (a,x). 
The most usual case that arises in quadrature applications is the case where w(x) = 1. In 
this case, we write 


Gif,ax, nt 1] =)) crea [soya 


i=0 


eit 
xX — Xj 
Ci = | dx 
7 Xj — Xj 
a j=0 J 


j#i 


x n a 
[ peoe-ctprexns =e ‘([]e-9] dx 


a a 
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and the nodes x; are the roots of the Legendre polynomial Py, (x) on [a, x] . Gaussian quadra- 
ture with the unit weight function is known as Gauss-Legendre (GL) quadrature. 


For the sake of clarity in the next section, we will replace the integration variable x with ft, 
while the upper bound in the integral remains as x. Also, we will denote the nodes by tj. 
So, we have 


: (2n+2) n(n : 
[ fod Gtf..axn03)-2 { on) dt (1) 


2.1.1 Gaussian Quadrature in terms of a Stepsize 

From this point onwards, we consider the specific case of GL quadrature. The roots of 
Pri (x) on [a,x] are symmetrically distributed within [a,x], but they are not equispaced 
and they are not located at a or x (so GL quadrature is open). Nevertheless, it is possible to 
write GL quadrature in terms of a stepsize, as can be done for Newton-Cotes quadrature. 


We define 
x-a 


n+2° 


so that h is the average separation of the nodes on [a, x]. We may now write 


G[f,1,a,x,n+1] = hy») dif (xi), 
i=0 


where 
(n+ 2) ¢; 


x-a 


dj = 


Furthermore, with the substitution t = a+ sh, where s € [0,n+ 2] is a continuous variable, 
and ¢; = a+ ojh, where 0; is an appropriate constant, we find 


so that 


x n+3 ¢(2n+2) n+2 a 2 
[f0e-ctirasney=- LEE foo} a @ 
a 0 


from which we conclude that GL quadrature is of order 2n + 3. 


The nodes and weights for GL quadrature for the interval of integration [-1, 1] are well- 


documented; a useful source is [5] (the Legendre polynomials can be generated using Bon- 
net’s recursion given in the Appendix. Thereafter, the roots can be found using numerical 
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software). It is easy to transform them to the case of an arbitrary interval [a, x]. If t; is a 
node on [-1, 1], then the corresponding node on [a, x] is given by 


t; = (i+1)(—*) +0 
and the relationship between the corresponding weights is given by 


xXx-a\ _ 
c= ( 2 a 


where ¢; is the weight associated with ¢; on [-1,1]. Note that t; and cj have now been 
expressed linearly in terms of x. 


2.2 Runge-Kutta methods 


We will assume that the reader has sufficient familiarity with Runge-Kutta (RK) methods for 
solving ODEs numerically. The RK method that we will employ here is an explicit seventh- 
order method, which we denote RK7. We refer the reader to Butcher [6] and Hairer et al [7] 
for information regarding these methods. 


3 Theory 


We now have 
n 


cusraneene S(Ahar(Bea( 


i=0 2 


=a) 


which we shall denote by Gf (a, x, n) from now on. 


The integral in (2) will be denoted by 


2 
Y(n) = aa | (Ie 2) ds 


and is easily determined for any given value of n (see item #1 in the Appendix). 


So, (2) becomes 
fl f (t) dt - Gp (a,x, n) = WPM f2™) (7 (x) V(n). 


Differentiating with respect to x gives 


dG 
f (x tt —~ = en (0) 


dh2"*3 Jud gp?) (CZ) dc 
Pa di =) (7), 
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which can be rearranged to yield 


i 
fly- ZL 4 243 
dc _ YQ) ee : (0) a 


dx panes LOA) 


Ae es (n) f2"*2) ( (C) ( ene 
Y (n) — (2) , 


Solving this ODE on [0, x] with a suitable initial value gives C(x), which then allows the 
error term to be determined. 


(3) 


4 Implementation 


Implementing a RK method to solve (3) requires addressing several issues. 


1. The factor f"*9) (2) in the denominator could be zero, or uncomfortably close to zero. We 
solve this problem by introducing a monomial offset, as we suggested in [1]. Essentially, 


we define a new function 
Dx23 


g@)=fG)+ Gap 


where 
D=1+max tiie (x), 
[0,x] 


and we apply our algorithm to g (x), instead of f (x). This gives 


atx) (2n+2) dpen+3 
alg Y(n) —§ (Ce) dx 
dx ents dgn+2)(Cx) 
dy 


g(x) - FE-¥(n) gen) (C,) (2ERE*) 

TOR EIG) = 
where C, and Gy are quantities relevant to the use of g(x) in (2). The result of this 
computation can be manipulated to yield the error term corresponding to f (x). See item 
#2 in the Appendix for detail. 

2. The factor h2"*3 in the denominator could also be problematic for values of x close to a. 
To counter this, we choose an initial value xp reasonably far from a, and then apply the RK 
method from xp up to a chosen upper limit xy. This gives the error term on the interval 
[x0, xn] , but not on the interval [a, x9]. To find the error term on [a, xo] , we apply the RK 
method from x9 down to a (a “backward” computation, so to speak) This action requires 
the use of a limit b in the integral (where x < b < xy), rather than a, and a suitable 
initial value, which is not the same as the initial value used for the forward computation 
(we discuss these initial values in the next point). We provide some detail regarding this 
procedure in items #3 and #4 in the Appendix. The nett result is a piecewise construction 
of the error term over the intervals [a, x9] and [x9, xn]. 
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3. It is necessary to determine an initial value in order to implement the RK method. Sub- 
stituting xp into (1) gives 


- ae fem?) (Z (xp)) fon a 
fre dt - Gy (a, Xo; n) = oa i (t = | dt. (5) 


If we evaluate the integral on the LHS to a suitable accuracy, (5) becomes a nonlinear 
equation that can be solved (numerically, if necessary) for ¢ (xo). Of course, the nodes t; 
are the relevant GL nodes on [a, xo] . This gives the initial value (xo, ¢ (xo)) for the forward 
RK computation on [x9, xy]. A suitable initial value for the backward RK computation 
on [a, xo] can be found by changing the limits of integration in (5) to give 


b b 2 
_ (2n+2) us 
[ro dt — Gr (xo, b, n) = eee iG | dt. 


Here, the upper limit in the integral is now b, the lower limit is xp, and the nodes ft; are 
the relevant GL nodes on [x9, b]. The initial value (xo, C(xo)) obtained here is not the 
same as in the previous case (since the integration limits are different), and it is probably 
a good idea to distinguish them through the notation 


(x0, Ga(Xo)) and (x, Cp (Xo). 


Of course, even if it is necessary to work with g(x), instead of f (x), the same process 
can be applied. 


5 Examples 


By way of numerical examples, we consider 
x 
‘| t’ sin tdt = 2x sinx — (x* - 2) cosx -2 


0 
x 


j[efdaens 


0 


i 
jf sintat = 1-cosx 
0 


0<x<10 (> a=0,xy = 10). 


Higher derivatives of f (x) = x? sin x and f (x) = sin x have zeroes on [0, 10] so that the use 
of a monomial offset is justified for these cases. Also, we will use xp = 0.5 and b = 1 for all 
three examples. For the sake of economy of presentation, we restrict our calculations to GL 
quadrature with n = 2, although our algorithm is general with regard to n. 
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If 5¢ (x) denotes the global error in C(x) due to the RK computation, we have 


pene O22) (x) + 5 (x) Va) = BMY (n) (OM (Z )) + 8c (x) FO" (C(x) +.) 
= 3 (n N auimes (4 Z (x)) 
Fe 5¢ (x) RV (n 5 Awan (5 7 (x)) diees 


where we see that 


pry Cy Acai, (C (x)) 
serves as a conditioning function, possibly amplifying |5< (x)]. 


We define E;", as the maximum magnitude of the difference between the true GL error 
and the GL error computed via C(x) on [xo, xn], E”, as the maximum magnitude of the 
difference between the true GL error and the GL error computed via C(x) on [a, xo], and 


Ciiy as the maximum magnitude of h’"'*Y (n) f (2n+3) (7 (x)) on [xo, xn] - 


In Table 1, we show values for various parameters of significance for each example. For D 


and Cy’, we have rounded up to the nearest integer. We did not need a monomial offset for 
Tis. 
Table 1: Values of various parameters of significance, with n = 2. 
f (x) Ca (x0) Co (x0) D ENN Eno CoN 
x? sin x 0.2498 0.7493 112 36x10" 37x10“ 907 
ex 0.2519 0.7519 - 22e10 52K 10 1534 
sin x 0.2505 0.7510 2 53x10" 26x10" 10 


Figures 1 — 3 (see end of paper) pertain to the example f (x) = x’ sin x. In Figure 1, we show 
C(x) on [0, xo] and [xo, 10] . The discontinuity is simply due to the two different initial values 
Ca (xo) and Cy (xo). In Figure 2 we show 


10 
[fea- G[f, 1, 0, 10, 3] 
0 


and 
(6) 2 - 
[x t) dt -| G[f,1,0, 10,3] + sau al netic f([] 0 ds). 
i=0 


The first of these is the uncorrected GL error, and the second includes the error computed 
via C(x) as a correction term. There is a clear improvement in accuracy. 


We observe, however, that both the uncorrected and corrected errors generally increase 
with x in Figure 2. We believe that there are two reasons for this: (i) the global RK error 
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in C(x) is likely to increase with x, and (ii) the condition function shows a clear increase 
in magnitude with x (see Figure 2), which only serves to amplify the effect of 5; (x) on the 
overall result. Indeed, the corrected error in Figure 3 is increased by some three orders of 
magnitude near the end of the interval due to the value of the condition function. 


6 Conclusion 


We have determined the error in Gauss-Legendre quadrature by solving a suitable ODE for 
C(x) - the central parameter in Gauss-Legendre quadrature - using an explicit Runge-Kutta 
method. Our algorithm is robust in that it caters for singularities and near-singularities in 
the ODE. Numerical experiments have shown that correcting the Gauss-Legendre approx- 
imation using the error determined in this way improves the accuracy of the quadrature 
by as much as 12 orders of magnitude. Our study has exposed the existence of a condi- 
tioning function, which can amplify the Runge-Kutta global error in C(x). This, in turn, 
suggests the need to control the Runge-Kutta global error in a step-by-step manner, as we 
have investigated previously [8][9][10]. We reserve this line of inquiry for future research. 


We have only considered Gauss-Legendre quadrature here, wherein the weight function 
w (x) is unity. Using other weight functions may be as simple as replacing f (x) with w (x) f (x) 
as the function to be integrated. This may necessitate a change of variable, as in the case 
of the Chebyshev weight function, but this, too, will be a topic of future research. 


Appendix 


1. The integral Y (n) is easily determined. We demonstrate for Y (2). With n = 2 we have 
n+1= 3nodes. These are the roots of the Legendre polynomial P3 (x) on [-1, 1], and are 


known to be 
- ne 
ty € 4-4/-, 0, Saat a 
5 5 


Using t; = -1+ ojh and h = 1/2 we find 


Hence, 


1 


‘ ae 2 
Y (2) = +I | a) ds 
= oy 
= 0.0081269841269804750372163226757039, 


and similarly for other values of n. As a matter of interest, we also find 


Y (3) = 0.0010984263083575790571899677416923. 
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2. The risk of a singularity in (3), due to the possibility of f?"**) (x) = 0 in the denominator 
on the RHS, is easily mitigated. Since we know f("*9) (x), it is easy to find a constant 
D such that f?"*9) (x) + D is not close to zero anywhere on [a, x] for all values of x of 
interest. Indeed, a suitable choice for D is 


D=1+max ae (x). 
[0,x] 


Define oe 
Dxe" 
B®) =fO)+ oy 
Hence, 
x x x on+3 
Dt*"™ 
t) dt = t) dt —— dt 
fea fpmas fm 
a a a 
and GL applied to g(x) is 
= 2n+3 
Gg (a,x,n) = G iz Cua | 
= 2nt3 
= Gr(a,x,n) + G|~ 1, a,x, n4 1). 
(2n + 3)! 
This gives 


fs@a-Gaxm= f pord-Glaxm 


a 


x 
Dte"*3 Dx23 


dt 
(2n + 3)! (2n + 3)! 


sane), 


so that 


x 


[Qa Elaxm- [sat Glaxm 


a 


x 
- | att 6] Pans 1 
(2n + 3)! (2n + 3)! 
a 
The first term in parentheses on the RHS is the error term that we obtain when we 
apply our algorithm to i * g(t) dt, and the second term in parentheses is easy to compute 
exactly. The net result is the error term corresponding to f (x), as desired. A similar 
process is used on the interval [x, b] . 
3. An easy way to implement an RK method 


Vit = Vi + hiss F (xi, Vir his1) 
hina = Xin — Xj (6) 
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in the negative x direction, starting at a point xo, is to simply label the nodes according 
to 


{. + 5X4, X3, X2, X14, Xo}, 


where we have xp > x1 > X2 > %3 > x4 > .... Note that, in (6), his; < 0 since xj41 < Xj. 
There are two points that must be considered regarding the backward RK computation: 
(i) we must use Gy (x, b, n) or Gy (x, b, n) for the GL quadrature term and, (ii), we must 
use —f (x) and —g (x) in (3) and (4), respectively, due to x being the lower limit. 

4. Assume a < z < x. We use the symbol Ar (x, 8, n) to denote the error in the GL quadra- 
ture Gf (x, B, n). We have 


froa-froas fro 


and 
b 
J Petrae = Gp(a.b,0) + Bp(a.b,1) 
fro dt = Gf (a, z,n) + Ar (a, z, n) 
: b 
[ fde= Gp G.b.m) Bp b0). 
Hence, 


Ag (a, z, n) = Gy (a, b, n) + Ag (a, b, n) - Ge (a, z, n) — Ge (z, b, n) - Ap (z, bn) 


which is easily determined since Ar (a, b, n) has been found via the forward RK compu- 
tation, Ar (z, b, n) has been found via the backward RK computation, and the other three 
terms are quadrature terms determined by direct substitution of relevant known values. 
This process can also be applied when g(x), instead of f(x), has been used - simply 
replace the subscript f with g in the above - to obtain A, (a, z, n). 

5. Bonnet’s recursive formula for the Legendre polynomials on [-1, 1] is 


7 (2n + 1) xP, (x) -— nPp-1 (x) 


Prit (x) n+i 


with Po (x) = 1 and P, (x) = x. 
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Figures 


For the convenience of the reader, all figures have been collected here. 


Forward RK 


—— Backward RK 


o) 


10 


Figure 1. C(x) for the forward and backward RK computations, 
for f (x) = x” sinx. The point xo is indicated. 


uncorrected 
04 corrected 


GL quadrature error 
be 
f 


Figure 2. Quadrature error (blue) and corrected quadrature 
error (green), for f (x) = x? sin x. Vertical axis is base-10 loga- 
rithmic. 
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ales) 


3 + | —— Condition function 


Condition function 


Figure 3. Condition function on [xo, xy] for f (x) = x? sin x. Ver- 
tical axis is base-10 logarithmic. The condition function clearly 
increases in magnitude quite significantly with x. 


